1. Field of the Invention
The invention relates generally to a cone beam computed tomography (CT) imaging system, and more specifically to a method and apparatus for identifying and correcting image inaccuracies caused by simplified processing of masked cone beam projection data.
2. Description of the Prior Art
Recently a system employing cone beam geometry has been developed for three-dimensional (3D) computed tomography (CT) imaging that includes a cone beam x-ray source and a 2D area detector. An object to be imaged is scanned, preferably over a 360.degree. angular range and along its entire length, by any one of various methods wherein the position of the area detector is fixed relative to the source, and relative rotational and translational movement between the source and object provides the scanning (irradiation of the object by radiation energy). The cone beam approach for 3D CT has the potential to achieve 3D imaging in both medical and industrial applications with improved speed, as well as improved dose utilization when compared with conventional 3D CT apparatus (i.e., a stack of slices approach obtained using parallel or fan beam x-rays).
As a result of the relative movement of the cone beam source to a plurality of source positions (i.e., "views") along the scan path, the detector acquires a corresponding plurality of sequential sets of cone beam projection data (also referred to herein as cone beam data or projection data), each set of cone beam data being representative of x-ray attenuation caused by the object at a respective one of the source positions.
As well known, and fully described for example in the present inventor's U.S. Pat. No. 5,257,183 entitled METHOD AND APPARATUS FOR CONVERTING CONE BEAM X-RAY PROJECTION DATA TO PLANAR INTEGRAL AND RECONSTRUCTING A THREE-DIMENSIONAL COMPUTERIZED TOMOGRAPHY (CT) IMAGE OF AN OBJECT issued Oct. 26, 1993, incorporated herein by reference, image reconstruction processing generally begins by calculating Radon derivative data from the acquired cone beam data. The Radon derivative data is typically determined by calculating line integral derivatives for a plurality of line segments L drawn in the acquired cone beam data. In the embodiment described in detail in the U.S. Pat. No. 5,257,183, Radon space driven conversion of the derivative data is used to develop an exact image reconstruction of a region-of-interest (ROI) in the object. Calculation of the line integral derivative data is a relatively complex and computationally time consuming task.
A cone beam data masking technique is known which improves the efficiency of the calculation of the derivative data in such Radon space driven reconstruction processing, as described in the present inventor's U.S. Pat. No. 5,504,792 entitled METHOD AND SYSTEM FOR MASKING CONE BEAM PROJECTION DATA GENERATED FROM EITHER A REGION OF INTEREST HELICAL SCAN OR A HELICAL SCAN, issued Apr. 2, 1996, also incorporated herein by reference. The masking technique facilitates efficient 3D CT imaging when only the ROI in the object is to be imaged, as is normally the case. In the preferred embodiment described therein, a scanning trajectory is provided about the object, the trajectory including first and second scanning circles positioned proximate the top and bottom edges, respectively, of the ROI, and a spiral scanning path connected therebetween. The scanning trajectory is then sampled at a plurality of source positions where cone beam energy is emitted toward the ROI. After passing through the ROI, the residual energy at each of the source positions is acquired on an area detector as a given one of a plurality of sets of cone beam data. Each set of the cone beam data is then masked so as to remove a portion of the cone beam data that is outside a given sub-section of a projection of the ROI in the object and to retain cone beam projection data that is within the given sub-section. The shape of each mask for a given set of cone beam data is determined by a projection onto the detector of the scan path which is above and below the source position which acquired the given set of cone beam data. The masked (i.e., retained) cone beam data is then processed so as to develop line integral derivative reconstruction data. An exact image of the ROI is developed by combining the reconstruction data from the various source positions which intersect a common integration plane. Hence, the masks are commonly referred to as "data-combination" masks.
Although the use of the data combination masks significantly reduces the processing required in Radon driven approaches, calculation of the line integral derivative data is still a relatively complex task and computationally time consuming. One known technique for developing the line integral derivative reconstruction data in such a Radon space driven reconstruction processing approach, is to use linograms. Although the linogram technique provides for a much faster and more simplified processing of the masked data for developing the line integral derivative reconstruction data, its creates less than exact, i.e., quasi-exact, reconstructed images.
Data-combination masks can also be used to improve the efficiency of the calculation of the derivative data in detector data driven techniques, such as those using Filtered BackProjection (FBP) techniques. A "simplified" FBP technique is described in the present inventor's U.S. Pat. No. 5,881,123 entitled SIMPLIFIED CONE BEAM IMAGE RECONSTRUCTION USING 3D BACKPROJECTION, issued Mar. 9, 1999, also incorporated herein by reference. This simplified technique reconstructs the image using 2D approximation data sets formed by ramp filtering of the masked cone beam data. The filtering is carried out in the direction of the projection of a line drawn tangent to the scan path at the source position that acquired that set of cone beam data. Although this technique is also less complex than the prior techniques, the reconstructed image is also quasi-exact.
Accordingly, the present inventor's U.S. Pat. No. 6,084,937 entitled ADAPTIVE MASK BOUNDARY CORRECTION IN A CONE BEAM IMAGING SYSTEM, issued Jul. 4, 2000, and also incorporated herein by reference, describes a technique for computing 2D correction data which, when combined with the ramp filtered 2D approximation data sets, is intended to yield an exact image reconstruction. The 2D correction data basically comprises a point spread function representative of image reconstruction processing for each point on the detector which intersects the boundary of the data-combination mask.
Although this technique, as well as the technique of the forenoted U.S. Pat. No. 5,504,792 are intended to yield and exact image reconstruction, the present inventor has realized that such techniques are in fact also quasi-exact. More specifically, in an exemplary filtered backprojection (FBP) image reconstruction, on a detector let u and v be the Cartesian coordinate axes of the detector with the v axis coinciding with the rotation axis, and let L(.theta.,s) denote the line on which projection/backprojection is carried out, where (.theta.-.pi./2) is the angle line L(.theta.,s) makes with the u axis, and s is the displacement of the line from the origin. Filtering of the cone beam image consists, either explicitly or implicitly, of the combined operation D.sub.t H, where D.sub.t is the differentiation (spatial) operation in the projected scan path direction t, and H is a shorthand notation for .intg.B(.theta.)D.sub.s (.theta.)P(.theta.)d.theta.. P(.theta.) is the 2D projection operation (line integration) at angle .theta., D.sub.s (.theta.) is the 1D differentiation operation with respect to s for the projection at angle .theta., and B(.theta.) the 2D backprojection operation at angle .theta..
FIGS. 1A and 1B illustrate this combined operation in an FBP reconstruction processing technique, such as described in the forenoted U.S. Pat. No. 6,084,937. FBP image reconstruction consists of two different kinds of processing: the first kind is 2-dimensional (2D) and the second kind is 3-dimensional (3D). In the 2D step each cone beam projection image is processed in a 2D space, and in the 3D step the processed 2D images are backprojected into a 3D object space. The 2D step consists of the following 4 FBP image reconstruction sub-steps for processing the cone beam data acquired at each of a plurality of the source positions (S.sub.i) along the scan path:
1. Compute a 1D projection (i.e., a line integral) of the cone beam image acquired on a detector plane 100, at each of a plurality of angles .theta.. This step is illustrated in FIG. 1A for a given angle .theta..sub.1 of a plurality of angles .theta.. A 1D projection 102 is shown at coordinates s, .theta..sub.1 comprising the integrated values of the cone beam image 104 on detector plane 100 along a plurality of parallel lines L(s, .theta.) that are normal to angle .theta..sub.1, each line L being at an incremental distance s from an origin O. As shown and described below in conjunction with FIGS. 2 and 3 illustrating the concept and implementation of data combination, the lengths of the lines L will be limited using the forenoted masking techniques. Generally, if the detector plane 100 comprises an N by N array of pixels, then the number of angles .theta. is typically given by .pi.N/2. PA0 2. Filter (differentiate) each 1D projection 102 in accordance with a d/ds filter, resulting in a new set of values at each of the s, .theta. coordinates, such as shown by the derivative projection 106 for the angle .theta..sub.1 shown in FIG. 1A. Note, the sum (integration) of the resulting values at these s, .theta. coordinates yield a quantity proportional to the Radon derivative for an integration plane Q(s, .theta.), as described above for Radon space driven image reconstruction processing. PA0 3. As illustrated by FIG. 1B, backproject the derivative projection 106 from each angle .theta. into a 2D object space 107 (which coincides with the detector plane 100). Lines 108 are representative of this backprojection, and spread the value from each s coordinate into 2D space 107 in a direction normal to each .theta., thereby developing contributions to a backprojection image 109. Note, 2D object space 107 has a size corresponding to a virtual detector which is enlarged (compared with detector having a height corresponding to the pitch of the scan path), so as to cover an entire region of interest (ROI in the object. This enlargement is required because the calculated Radon data affects the reconstruction of the entire Q plane, and not just the partial plane represented by the data combination mask. PA0 4. Perform a 1D d/dt filtering of the backprojection image formed in 2D space 107 by step 3. The 1D filtering is performed in the direction of the scan path, i.e., along lines 110, where t points in the direction of the projection of a line drawn tangent to the scan path.
In a further "simplified" FBP reconstruction processing technique, the above steps 1-3 are combined into a single step of "Feldkamp" ramp filtering of the masked 2D projection data in the t (projection of a tangent to the scan path) direction, as described in detail in the present inventor's forenoted U.S. Pat. No. 5,881,123.
After step 4 above, and as shown in FIG. 1B, the 3D step comprises performing a weighted 3D backprojection of the thus filtered data from 2D space 107 (i.e., from each pixel in the detector) onto a plurality of sample points P in a 3D object volume 112. The density assigned to each point P is weighted by the inverse of the square of the distance between the sample point and the location of the x-ray.
Common in image reconstruction techniques using a data combination mask, such as in the above described FBP and simplified FBP techniques, is the forenoted data combination masking process for limiting the x-ray coverage of the integration plane at each source position to the angular range bounded by the prior source position below and the subsequent source position above the current source position. Such data combination from several source positions which intersect a common integration (Q) plane is illustrated in FIG. 2. More specifically, FIG. 2 represents a typical integration plane Q(s, .theta.) intersecting a cylindrical object and a spiral scan path, which is assumed to wrap around the object on an imaginary cylinder and having top and bottom circular scan paths even with top and bottom edges of a region-of-interest (ROI) in the object. An edge view of plane Q is illustrated in FIG. 4 (described in the detailed portion of the description). Since a non-vertical plane will intersect a cylinder in an ellipse, the plane Q(s, .theta.) intersects the object and the cylindrical spiral scan path in two ellipses, E1 and E2, respectively, one inside the other, as shown in FIG. 2.
Since the spiral path lies on the scan path cylinder, it intersects the plane Q in points that lie on the ellipse E2. These source positions are shown as S1, S2, and S3 in FIG. 2. Similarly, it is easy to see that the top scan path circle intersects the plane in two points T1 and T2 which lie at the intersection between E2 and the top edge of the object's ROI and that the bottom circle intersects the plane in the two points B1 and B2 which lie at the intersection between E2 and the bottom edge of the object's ROI. Other integration planes may have more or less spiral scan path intersections, depending upon their orientation, and may not intersect either of the top or the bottom circle scan paths.
As is apparent from FIG. 2, the source positions which illuminate that portion of integration plane Q that lies within the ROI are T.sub.2, S.sub.1, S.sub.2, S.sub.3, and B.sub.2. Complete X-ray coverage of region-of-interest 200 of this portion of the integration plane can be achieved by suitably combining the data acquired at these 5 source positions, as indicated in FIG. 3. For example, at T.sub.2 we only use the cone beam data within the angle bound by T.sub.1 T.sub.2 and T.sub.2 S.sub.1, and at S.sub.1 we only use the cone beam data within the angle bound by T.sub.2 S.sub.1 and S.sub.1 S.sub.2. And so on. Five partial planes P1 through P5 are therefore defined by the source positions T.sub.2, S.sub.1, S.sub.2, S.sub.3, and B.sub.2, which do not overlap and together completely cover the portion of plane Q that lies within the ROI. In this way the totality of the cone beam data from each of the contributing source positions illuminates the entire Q plane only once without any overlap. Further details of this data combination technique can be found in the present inventor's earlier cone beam patents, such as U.S. Pat. No. 5,463,666 or U.S. Pat. No. 6,084,937).
The mask consists of a top curve and a bottom curve formed by projecting on to the detector plane the spiral turn above and the spiral turn below the current source position. For a flat detector located at the rotation axis such that the line connecting the source to the detector origin is normal to the detector plane, the equation for the top curve for the spiral scan is given by: ##EQU1##
where R is the radius of the spiral and h is the distance between adjacent spiral turns (the pitch). The bottom curve is the reflection of the top curve about the origin, i.e. (u,v) {character pullout} (-u,-v). The shape of one such spiral mask is shown in FIG. 3. Further detail about the generation and use of data combination masks can be found in my prior cone beam patents, such as U.S. Pat. No. 5,463,666 or U.S. Pat. No. 6,084,937.
The present inventor has realized that although his forenoted U.S. Pat. No. 6,084,937 entitled ADAPTIVE MASK BOUNDARY CORRECTION IN A CONE BEAM IMAGING SYSTEM is intended to reconstruct exact images, in fact there are still image errors. This is because the masking process is only an approximation for limiting calculation of the line integrals intersecting the cone beam image in order to achieve proper data combination. Thus, the reconstruction processing is in fact quasi-exact, i.e., the reconstructed images have artifacts.
It would be desirable to provide a way to easily identify and prevent/correct such image artifacts when performing a simplified image reconstruction processing.
The present inventor has realized one source of image errors comes from "second intersections" between the integration line segments and the mask boundary. More specifically, consider the top mask boundary and the line L illustrated in FIG. 3, where the spiral path which projects onto the mask boundary scans from bottom to top in a clockwise direction. Line L intersects the top mask boundary at 2 points M.sub.1 and M.sub.2. It then follows that the spiral scan path intersects the integration plane Q defined by the line L and the current source position in the following order: current source position, followed by M.sub.1, then followed by M.sub.2. Thus the portion of the line that conforms to data combination, i.e. the x-ray data in the angular range bounded by the previous source position below and the subsequent source position above, is the segment to the right of M.sub.1. When filtering the cone beam image, it is this segment alone that should be included in the projection operation P(.theta.), i.e. the line integration operation. However, in the above-mentioned filtered-backprojection reconstruction algorithms (as well as in the above noted linogram technique known for use in Radon space driven reconstruction processing) it is assumed that the entire portion of the line inside the mask is included in the projection operation. That is to say not only the segment to the right of M.sub.1 but also the segment to the left of M.sub.2 is included in the projection operation. Such unneeded contribution to projection occurs wherever: (1) the line of integration intersects the same mask boundary twice, and (2) the second intersection point lies within the detector. More specifically, for the top and bottom mask boundary, errors are caused by line integration on the line segment to the left and right, respectively, of the second intersection point. It is noted that in the exact Radon space driven reconstruction techniques noted above, these second-intersection contributions were avoided via explicit tabulation for each line integral. Such tabulation is not practical for simplified Radon driven processing (i.e., the linogram technique or for filtered-backprojection processing).
It would be desirable to provide a way to easily identify and prevent/correct such image artifacts due to second intersections with the data combination mask when performing simplified image reconstruction processing.